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Abstract 

We propose two different macroscopic dynamics to describe the decay of 
metastable phases in many-particle systems with local interactions. These 
dynamics depend on the macroscopic order parameter m through the re- 
stricted free energy F(m) and are designed to give the correct equilibrium 
distribution for m. The connection between macroscopic dynamics and the 
underlying microscopic dynamic are considered in the context of a projection- 
operator formalism. Application to the square-lattice nearest-neighbor Ising 
ferromagnet gives good agreement with droplet theory and Monte Carlo sim- 
ulations of the underlying microscopic dynamic. This includes quantitative 
agreement for the exponential dependence of the lifetime (r) on the inverse 
of the applied field H, and the observation of distinct field regions in which 
A = dln(T) / d\H\ l ~ d depends differently on \H\. In addition, at very low tem- 
peratures we observe oscillatory behavior of A with respect to \H\, due to the 
discreteness of the lattice and in agreement with rigorous results. Similarities 
and differences between this work and earlier works on finite Ising models in 



the fixed-magnetization ensemble are discussed. 
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I. INTRODUCTION 



Metastable phases are observed in a wide variety of systems that exhibit first-order 
phase transitions. A few examples are supercooled fluids, permanent magnets, ferroelectrics, 
certain alloys, the "false vacuum" associated with the electroweak phase transition, and the 
supercooled quark/gluon plasma associated with the QCD confinement transition. In recent 
decades, much attention has been focused on the study of metastable phases and the rate 
at which they decay to thermodynamic equilibrium, but a fully satisfactory description 
has remained elusive. A recent review with numerous references to specific realizations of 
metastable behavior in real and model systems is found in Ref. [l]. 

In certain systems with weak long-range interactions, infinitely long-lived metastable 
phases can exist in the thermodynamic limit [Q. However, in systems with short-range 
interactions, there exist no such stable non-equilibrium states, even in the thermodynamic 
limit. Nevertheless, for large but finite systems, the relaxation time for short-range models 
can be extremely long compared with any finite observation time 0-0. Here, we define the 
term "metastability" to include this phenomenon in short-range models. The long relaxation 
time is mainly due to the large free energy of the local fluctuations that must spontaneously 
arise in order for the system to decay into a globally stable phase. Due to the long relaxation 
time, it is difficult to tell metastable phases from globally stable ones by observing only 
short-time fluctuations. The explorations of phase space, characteristic of the metastable 
phase, are expected to be those included in a constrained partition function that excludes 
the microstates that dominate in equilibrium PJ7|]. The application of such ideas to a field- 
theoretical droplet model with Fokker-Planck dynamics has shown that close to coexistence, 
the nucleation rate for droplets of the equilibrium phase is proportional to the imaginary 
part of a complex-valued constrained free energy obtained by analytic continuation from the 
equilibrium phase into the metastable phase [jSHTOfl. 

Recently, complex-valued constrained free energies were numerically obtained for both 
the two-dimensional nearest-neighbor Ising ferromagnet fn|,[T2| and for models with weak 
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long-range forces [p~3| pT5f| by a constrained-transfer-matrix method introduced by one of us 
fl6|| . Although no dynamical aspects were explicitly considered to obtain the constrained free 
energies, the average free-energy cost of a critical droplet was obtained over a wide range of 
fields and temperatures, in good agreement with the predictions of field-theoretical droplet 
models [^-[TD|,|T^,[T3|,|T^| and Monte Carlo (MC) simulations These results indicate 



the relevance of purely static properties, such as the free energy, to the relaxation behavior 
of metastable phases. Whereas any physical dynamic (consistent with a real experimental 
situation) is bound to give the correct equilibrium Boltzmann distribution for an infinitely 
long observation time, it is not yet clear how relevant the static properties of a model are 
to the study of the dynamical relaxation of a metastable phase towards equilibrium. 

The observations discussed in the preceding paragraph raise the interesting possibility 
that the information stored in static quantities may be sufficient to describe the salient 
features of the relaxation behavior of a metastable phase, even in a short-range-force system. 
A quantity which contains all thermodynamically relevant equilibrium information is the 
restricted bulk free energy, 

F(m) = F (m) - (3HNm , (1) 

where m is the macroscopic order parameter conjugate to the external field H. Here (3=1 /T 
is the inverse temperature with Boltzmann's constant ks—l, and (3 has been absorbed in 
F(m) and F Q (m). One can obtain F(m) either exactly from exact enumerations for small 
systems, or approximately, up to an additive constant, by Monte Carlo simulation. The 
importance of the detailed shape of the zero- field free energy F Q (m) for two-phase equilibria 
and nucleation barriers has previously been discussed by Binder and coworkers fi8H20|. 
Generalizations of Eq. ([I]) to consider several macroscopic densities and their conjugate 
fields are straightforward. 

Since the restricted bulk free energy is (by definition) projected onto a space spanned by 
one or a small number of macroscopic densities, all detailed information about microscopic 
spin configurations is lost. In this paper we investigate the possibility that macroscopic 



dynamical properties that are common to several different microscopic dynamics may never- 
theless be extracted from the information retained in F(m). For this purpose, we construct 
two different macroscopic, discrete-time dynamics, each defined by a separate master equa- 
tion for the order-parameter distribution function. Both master equations are subject to the 
following two restrictions. 

(1) The order parameter m is allowed to change only by a finite amount during each discrete 
time step. (Locality in m.) 

(2) The dynamics should reproduce the correct F{m) in equilibrium. (Correct static prop- 
erties of F(m).) 

Although the relevance to metastable decay of the Hohenberg-Halperin scheme of dynamic 
universality classes |21[ is not completely clear, we demonstrate in this work that the re- 



quirements (1) and (2) are sufficient to make our macroscopic dynamics consistent with 
microscopic dynamics in the class of Model A: systems with a nonconserved scalar order 
parameter and local dynamic. 

For short-range-force models, the sequence of microscopic configurations that constitutes 
a particular realization of a MC simulation cannot be deduced from the corresponding 
sequence of values of macroscopic variables. In contrast, for models in which each site 
interacts equally with all other sites while the total interaction energy remains independent 
of the system size (equivalent-neighbor models), all configurations with identical values of 
the order parameter are equivalent, so the dynamical properties can be exactly obtained 
from the restricted free energy, as has been shown by Griffiths, et al. [22|. In the equivalent- 
neighbor limit, one of the macroscopic dynamics that we propose in this work reduces to the 
Metropolis |23] version of the heat-bath dynamic studied in Ref. ||22|| . (For discussions of the 



distinctions between Metropolis and heat-bath or Glauber dynamics, see e.g. Refs. fl|,f24|j.) 
Since it is well known that both the equilibrium and the metastable properties of equivalent- 
neighbor models are exactly described by mean-field theory in the thermodynamic limit 
PJT3|-P^o]], these models are often referred to as "mean-field models" |22|| . Consistent with this 



usage, we call the class of dynamics that we define here "macroscopic mean-field dynamics." 



Our proposed dynamics may be considered as approximations to the dynamic one would 
obtain by projecting the microscopic dynamic onto a master equation for the macroscopic 
order-parameter distribution, using a projection-operator technique p5|-p0|]. This point of 
view is further explored in Appendix [A[ Since the decay of metastable phases is a nonlinear, 
nonequilibrium problem, it is worth considering the extent to which nonlinearities and corre- 
lations in the microscopic dynamic are included in the proposed macroscopic dynamics. We 
therefore point out that, although specifically nonequilibrium correlations are not included 
by virtue of the loss of spatial resolution resulting from the projection of the dynamic onto 
the macroscopic order parameter, equilibrium correlations are included through their effect 
on the highly nonlinear restricted free energy F(m). However, since in this paper we only 
consider a single macroscopic variable, effects of nonlinear interactions between macroscopic 
variables are not included. 



The rest of this paper is organized as follows. In Sec. [Q| we introduce our two spe- 
cific macroscopic mean-field dynamics, emphasizing their common physical motivation. In 
Sec. [TTI] we show how to obtain F{m) using MC simulations. In Sec. [TV] we summarize the 
relevant droplet-theory predictions for the decay of metastable phases. In Sec. |V| we apply 
our macroscopic dynamics to the relaxation of the metastable phase in the two-dimensional 
ferromagnetic nearest-neighbor Ising model below its critical temperature. In Sec. |Vl]we dis- 
cuss the connections between our results for the metastable lifetimes and the detailed shape 
of F(m), with particular reference to the earlier work by Binder and coworkers fi8H20|. 
Finally, in Sec. [VII] we summarize our results and discuss some implications of this study. 

II. MACROSCOPIC MEAN-FIELD DYNAMICS 

To set the stage for our study, we first consider the microscopic Metropolis dynamic for 
a mean-field Ising ferromagnet in which each spin interacts with equal strength with every 



other spin in the system. As was already shown by Griffiths et al. [22], this exactly defines 



a macroscopic dynamic which depends on the configurations of the system only through the 
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restricted free energy F(m). Next, for Ising models with finite interaction range, we propose 
two macroscopic dynamics, which can be constructed from F(m) for the corresponding 
model while satisfying the two conditions introduced in Sec. |: (1) locality in m and (2) 
correct static properties of F(m). In proposing these dynamics we make as few specific 
physical assumptions as possible beyond the conditions (1) and (2). 

The ferromagnetic Ising model with equivalent-neighbor interactions is defined by the 
Hamiltonian 

H=-(J/N)Y i 8 i 8 j -H'£ i 8 i , (2) 

i<j i 

where Sj=±l are N Ising spins, H is an external magnetic field, and the sums J2i<j an d J2i 
run over all N spins with l<i<j<N and l<i<N, respectively. For convenience we set the 
interaction constant J equal to 1. The magnetization per spin, 



m 



AT 1 ^, (3) 



is the order parameter conjugate to H, and the number of up spins is related to m as 

N 

n = — (1 +m) . (4) 



With these definitions, Eq. (g) can be written as |22 



U = E(n) = - {2n 2N N? - H(2n -N) + ±. (5) 

We consider a microscopic Metropolis dynamic in which the spin at a randomly selected 
site i is flipped from — > —s^ with probability 

p(x -> x') = exp [min {0, (3 (E(x) - E(x'))}} , (6) 

where E(x) and E(x') are the energies of the microscopic spin configurations x={sj} and x' 
before and after the flip, respectively. The microscopic detailed-balance condition is satisfied 
pi| . For this model, the one-step transition probabilities Wi(n,n') from states with order 



parameter n to states with n' are |2^ 



W 1 (n,n + 1)= (l-^)exp[jmn{0,P(E(n)-E(n+l))}] (7a) 

Tl 

Wi(n, n - 1) = — exp [min {0, (3 {E(n) - E{n - 1))}] (7b) 
W^n.n) = 1 - W 1 (n,n + 1) - W^n.n- 1) . (7c) 

Here the arguments of the matrix elements Wi(n, n') are all between and N, and ma- 
trix elements with arguments outside this range are identically zero. In this dynamic the 
probability for choosing an up (down) spin is n/N (1 — n/N). 

Since the value of n in the equivalent-neighbor model uniquely specifies the energy of the 
spin configuration, j3E(n) can be replaced by F(n) + S(n), where F(n) is the restricted free 
energy [F(n) = (3E(n) — S(n) with (3 absorbed in F], and S(n) = lnQ(n) is the Boltzmann 
entropy for the density of states Q(n) = N\/n\(N — n)\. The probability that the system 
has n up spins is proportional to exp[— F(n)] Eq. (^) thus becomes 

Wi(n, n + 1) = ^1 - exp [min {0, (F(n) - F(n + 1) + S(n) - S{n + 1))}] (8a) 

Tl 

W x (n,n-\) = —exp[min{0,{F(n)-F{n-l) + S{n) - S(n-1))}} (8b) 
Wxin, n) = \-W x [n,n+\)- W^n, n - 1) . (8c) 

Since exp[S'(n) — S(n — 1)] = (N — n + l)/n, it is straightforward to show that 

^| = exp[F W -F M ], (9) 

which can be considered to be the macroscopic detailed-balance condition between states 
with order parameter n and n'=n±l. Consequently, the equilibrium probability distribution 
for the order parameter n resulting from this Metropolis dynamic is correctly proportional 
to exp[— F(n)]. 

Once the transition probability matrix Wi(n,n'), which for iV Ising spins is an (N + 
1) x (N + 1) tridiagonal matrix, is constructed from F(n), particular realizations of this 
stochastic process can be created, starting from an arbitrary initial state. The average first- 
passage time (r) to the globally stable phase, starting from a metastable phase, can be 
easily calculated from W\{n, n') by the methodology of absorbing Markov chains |32| (see 
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Appendix |B|) . Higher moments of the first-passage time r can also be obtained by the same 
method. We again emphasize that for the equivalent-neighbor model, the time evolution of 
the macroscopic order-parameter distribution given by Eq. (|8|) is an exact consequence of 
the underlying microscopic dynamic p2| . 

Next, as a prototype ferromagnet with short-range interactions we investigate the 
nearest-neighbor square-lattice Ising model for which the Hamiltonian is 

U = -J E s l s J -HY,s l , (10) 

<i,j> i 

where Si = ±1 is the Ising spin at site i. The interaction constant J will be set to 1 as before, 
periodic boundary conditions are used, and H is the applied field [Q. The sums J2<ij 



> 



and J2i run over ai l nearest-neighbor pairs and over all N=L 2 sites on a square lattice. In 
contrast to infinite-range models, for models with finite interaction range, such as given 
by Eq. (|T0|), Eq. ([8]) can not be exactly derived from the microscopic Metropolis dynamic 
defined by Eq. (^), even if F(n) is numerically calculated to give the correct functional form 
appropriate to the particular model. The reason for this is that the local environment with 
which a particular spin interacts is no longer uniquely determined by the macroscopic order 
parameter. Nevertheless, we can define a macroscopic dynamic for a particular short-range 
model through Eq. (f|) with the appropriate form for F(n) and consider it as an approxima- 
tion for the true order-parameter dynamic observed in a microscopic MC simulation. This 
approximate macroscopic dynamic is our "mean-field dynamic No. 1" (MFD1). Although 
MFD1 is not exactly derived from any particular microscopic dynamic, it satisfies the two 
conditions of locality in the order parameter m and correct static properties of F(m), which 
we introduced in Sec. Q. 

We emphasize two important features of MFD1. First, regardless of the functional form 
of F{m) for the particular microscopic model, this form is correctly reproduced by the 
dynamic. Second, MFD1 is not the only macroscopic dynamic capable of correctly yielding 
F(m). This is analogous to the fact that there exist many different microscopic dynamics 
which can be successfully used to study the equilibrium properties of a single model. In fact, 
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the macroscopic detailed-balance condition, Eq. @, is sufficient to ensure that a dynamic 
yields exp[— F(m)} as its equilibrium distribution. Perhaps the simplest Metropolis- type 
dynamic which satisfies both Eq. (|9[) and our requirement of locality in m, is defined by the 
transition-matrix elements 

W 2 (n, n + 1) = - exp[min{0, F(n) - F(n + 1)}] (11a) 
W 2 (n, n - 1) = ^ exp[min{0, F(n) - F(n - 1)}] (lib) 
W 2 {n, n) = l- W 2 {n, n + 1) - W 2 {n, n - 1) , (11c) 

where the range of the arguments is the same as for W\, and W 2 =0 for out-of-range argu- 
ments. These transition probabilities define our "mean-field dynamic No. 2" (MFD2). 

In Sec. [V] we compare the macroscopic order-parameter dynamic observed in MC simu- 
lations of the two-dimensional Ising ferromagnet with the microscopic Metropolis dynamic 
defined by Eq. (|B|) to our approximate macroscopic dynamics, MFD1 and MFD2. This allows 
us to investigate the relevance of the equilibrium properties of a system to the dynamical 
relaxation behavior of its metastable phases. 

III. RESTRICTED FREE ENERGY BY MONTE CARLO SIMULATIONS 

In this section we describe how to obtain, from microscopic, equilibrium MC simula- 
tions, the bulk restricted free energy F(m), which is defined through the restricted partition 
function 

exp[— F{m)] = ^<5(m(x) — m) exp[— /3E(x)] , (12) 

X 

where the sum is over all possible microscopic spin configurations x, and 5 is the Kronecker 
delta function. It is straightforward to show that the probability distribution for m is 
proportional to exp[— F(m)] |3T|]. For Ising models with short-range interactions below 
the critical temperature T c , F (m) (i.e. F(m) for H=0, as defined in Eq. (|1|)) has two 
symmetrical minima, and the bulk free-energy barrier separating these minima diverges as 
the linear system size L for L^>1 in two dimensions (as L d ~ l in d dimensions) [j3lH . 
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Recently, significant progress has been achieved in the search for more efficient MC 
sampling algorithms for systems in which different subsets of phase space are separated by 



large free-energy barriers p4|-[41|j. Here we use a variation of the multicanonical method 



D), employing the notation of Ref. JTT|. 
The partition function for H=0 and inverse temperature (3 is 

Z(P,H=0) = E E , M ex P [S(E,M)-pE} = Em exp[-F (M)] 



(13) 



where 



exp[-F (M)] = ]T exp[S(E, M) - (3E] 



(14) 



exp[S(E, M)]=Q(E, M) is the density of states, and E and M=mL 2 are the bulk internal 
energy and magnetization, respectively 

The detailed-balance condition for the MC simulation can be written as 



W(x —> x') 



exp 



-(3{E(x) - E(x)} - [j(M(x')) -J(M(x))} 



(15) 



W(x' — > x) 

where J can be any arbitrary function of M. For an ergodic MC algorithm, the resulting 
distribution (histogram) of the sampling has been shown [il to be 



H(E, M) oc exp[S(E, M) — (3E - J{M)\ . 



(16) 



From Eqs. (|T4|) and (|T^) we get 



Therefore, 



H(M) = Y,H(E,M) oc exp[F (M) - J(M)\ 

E 



F (M) = J{M) + In H{M) 



(17) 



(18) 



up to an additive constant. 

As in Ref. ||T|, the quantity Fq(M) can be obtained using Eq. flT8] ) in an iterative 
fashion. The old estimates of Fq(M) and J(M) can be used in a MC procedure to obtain 
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a new histogram and consequently a new and better estimate for Fq(M) from Eq. (|18|). 
Although this approach is more efficient than conventional sampling methods ||40|| , we find 
that it does not suffice to obtain Fq(M) at very low temperatures. This is due to the large 
size of the exponent in Eq. (O) for large j3. 

As in Ref. |41|| , we instead sample the Boltzmann entropy S(E, M) directly by imposing 



the detailed-balance condition 

W(x — >■ x' 



W(x' — > x) 
which yields the histogram 



exp 



J{E(x),M(x)} - J{E(x'),M(x')} , (19) 



H(E, M) oc exp[S(£, M) - J(E, M)\ . (20) 
Therefore, S(E, M) can be obtained, up to an additive constant, by iteratively evaluating 

S(E,M) = J(E,M) +kiH(E,M) (21) 



where J(E, M) is an input and H(E } M) is the resulting histogram of J(E, M) from Eq. (|19|). 
Once S(E,M) has been obtained to the desired accuracy, Fq(M) can be obtained from 
Eq. (JTJ) for any (3. We were able to obtain F (M) for any (3>f3 c for the nearest-neighbor 
Ising ferromagnet on LxL square lattices with L<24. The if -dependent F(m) are trivially 
constructed from F (m) by Eq. ([]]). 

In Fig. we show F (m) for T=0.8T C w 1.8153 (Fig. 0(a)) and for T=0.2 w 0.0881T C 
(Fig. |](b)). The sawtooth-like behavior of Fo(m) for T=0.2 is due to the discreteness of 
the lattice. For T=0.8T C we used Eq. dl5|) to obtain F (m) for lattices up to L=64. The 
vertical arrows marked m c indicate the exactly known magnetization for which the most 
likely configuration with the given magnetization changes from a slab (for \m\ < m c ) to a 
droplet (for \m\ > m c ) in an infinite system |4"2"f . In Fig. |l](a) we have also marked the value 



of m c corresponding to L=24, taken from Fig. 3 of Ref. |ji2|| . For further discussion of the 
droplet-to-slab transition and its significance, see Sec. [VI]. 
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IV. DROPLET THEORY AND ITS PREDICTIONS 



In this section we review some predictions of droplet theory for the metastable lifetime 
(r) in (^-dimensional ferromagnetic systems with short-range interactions and local dynamics 
|33| . We leave details of the theory to Refs. and the references cited there. Starting 
with a large magnetization opposite to the applied field, we consider its relaxation and define 
(r) as the average time it takes for m(t) to reach a particular cutoff value, m cut . 

In addition to the lattice constant, which we take as unity, five length scales are important 
in the droplet theory for relaxation of metastable phases. These are the linear system size 
L, the radius of the critical droplet R c , the average distance between supercritical droplets 
Ro, and the single-phase correlation lengths in the stable and metastable phases, £ s and 
£ ms , respectively. Since the temperatures of interest in the present study are well below T c , 
£ s is practically field independent and of the order of unity, and since the field strengths 



considered are moderate, the same is the case for £ ms |12|]. Thus we are left to consider the 
interplay between three lengths: L, Ro, and R c , all of which are larger than unity in the 
temperature- and field regimes of interest. 

By comparing the bulk free energy gained by creating a droplet of the equilibrium phase 
in the metastable background with the free-energy cost of creating the droplet surface, one 
can show that the critical radius, beyond which the droplet is more likely to grow than to 
shrink, is 

= (d-l)go(T) ^ (d-l)tro(T) 

Am\H\ ~ 2m cq (T)\H\ 1 ' 

where Am is the magnetization difference between the metastable and the stable phase, and 
m eq (T) is the spontaneous equilibrium magnetization. Also, cr (T) is the equilibrium surface 
tension along a primitive lattice vector, and is assumed to be equal to the surface tension 
in the metastable phase. For the two-dimensional Ising model, both gq(T) and m eq (T) 



|U are exactly known. The relaxation proceeds in different ways, depending on the relative 



sizes of L, Ro, and R c . 
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If I if I is sufficiently small that R C >L, then the saddle-point configuration is a slab span- 
ning the system in d — 1 dimensions, and L is the most important length scale for the 
relaxation. The metastable lifetime is then determined by the surface free energy of such a 



slab. Consequently, for periodic boundary conditions it increases with L as 



(r(T,if,L))~exp(2/5(T (r)L rf - 1 ) . 



(23) 



This region of ultraweak fields is called the "coexistence" (CE) region ||, since the dynamic 
is similar to that at if =0, where two competing bulk phases coexist. 

As |if| is increased, R c becomes smaller than L. The crossover to the field regimes 
where the relaxation is dominated by critical droplets smaller than L has been called the 
"thermodynamic spinodal" (THSP) [§]]. The crossover field if thsp can be estimated by 
requiring that the critical droplet should occupy a volume fraction 



1 



m 



m eq (T) J 

corresponding to the cutoff magnetization m cu t. This yields pi 



(24) 



if 



THSP 



(0) = T 



(d-l)H(T) 



-i l/d 



(25) 



L [ 2m eq (T)0 _ 

where S(T) can be calculated from m eq and the anisotropic equilibrium surface tension by 
the equilibrium Wulff construction H43|j48| to obtain the critical droplet shape |TT| , p!2|j49[ . 

For fields somewhat stronger than if thsp 5 so that L 3>i? c 3>l, the nucleation rate per unit 
volume for critical droplets becomes p HlO| , p^ 



T(T,H) oc |ff| 6+c exp 



-/35(T) 



Iff 



^r{i + o(ff 2 )} 



(26) 



The exponent b is a universal exponent related to excitations on the droplet surface, and 
the nonuniversal exponent c gives the ff dependence of a "kinetic prefactor" wn i cn 
contains all dependence on the details of the dynamics. For d—2 and 3 it is expected that 
b=l and —7/3, respectively |T7|]. This has been confirmed by several methods, most recently 



by constrained-transfer-matrix calculations |[TT]JT^1 . For dynamics that can be described by 
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a Fokker-Planck equation it is expected that c=2 P, |10| , |T7|| . In a recent MC study for d=2 
it was confirmed that 6+c~3 for the Metropolis dynamic with updates at randomly chosen 
sites |J. 

If Rq^>L^$>R c , a single critical droplet is sufficient for macroscopic decay to occur before 
additional droplets nucleate. This region is called the "single-droplet" (SD) region [J5]. 
Assuming that the exact value of m cut is not an important factor (see the discussion in 
Ref. H), the average relaxation time (r) to m cut can be written as 



(r(T, H, L)) « [L d T(T, H)] oc L- d \H\^ b+c) exp 
This leads to 



§S{i+^ 2 )} 



A(T, H) = dH :^ n H : L)) = /35(T) + ^li^" 1 



(27) 



(28) 



where we have neglected higher-order correction terms. 

If L^>R ^>R C , many critical droplets nucleate before the decay of the order parameter 
can proceed to a macroscopic extent. This region is called the "multi-droplet" (MD) region 
0. In this region, we expect (r) to be independent of L. With the assumption that the 
radial growth velocity of the supercritical droplets is proportional to the applied field H 
)2|, (t) is predicted to be 1,|5|-|3 

(3E(T) 



, , , , , b+c + d 

{t(T,H)} ~ |#|"^+^ exp 



(d+l)\H\ 



-{l + 0(H 2 )} 



which leads to 



A(T,tf) = ^ + ^W- 



(29) 



(30) 



d+l ' d 2 -l 

where higher-order corrections again have been neglected. One also has the mean droplet 
distance p], ^B] , |^7 | 

(3E(T) 



Ro(T,H) oc \H\—*&-exp 



(31) 



(d+l)\H\ d - 

The field i^DSP, which separates the SD and MD regions was called the "dynamic spinodal 
field" in Refs. |H|. It can be estimated by setting Rq(T, H)(xL with a proportionality 
constant of order unity, which gives the scaling relation, 
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rr ( PZ(T) 

-"DSP 



1 

d-1 



(d+l)\nL / 



(32) 



InL J VlnL, 

expected to be asymptotically valid for nonzero T. 

These droplet-theoretical results can be summarized as follows. For a given system size 
L, if \H\ is small enough so that R C >L, then L is the most important length scale for the 
relaxation, and the system is in the CE region. As \H\ is increased beyond .HthsP) Rc 
becomes smaller than L. If Rq^>L^R c , the system is in the SD region, and the relaxation 
is characterized by Eq. (ffif), from which one can estimate b+c and H(T). As is increased 
still further, so that R c <^iRo<^L, the system is in the MD region. The relaxation is then 
characterized by Eq. (fjPD , from which one again can estimate b+c and S(T). The crossover 
field that separates the SD and MD regions, .£/dsp> a ls° separates the "stochastic region" 
and the "deterministic region" ||]. In the "deterministic" region, the average metastable 
lifetime (r) is quite short, and the standard deviation of r is much smaller than (r). In the 
"stochastic" region, the decay approximately follows a Poisson process, so that the standard 
deviation of r is comparable to (r). As \H\ becomes very large, the droplet picture becomes 
inappropriate, and Rq and R c become comparable to the lattice spacing. The lifetime is then 
on the order of one MC step per site (MCSS), and this region is called the "strong- field" 
(SF) region. The crossover between the MD and SF regions is marked by the "mean-field 



spinodal field" -£^mfsp H; which can be estimated as the field at which 2R C =1 |fL2] . In Fig. |2| 
we show a schematic diagram for the relaxation behavior, illustrating the four sub-regions 
of characteristic relaxation behavior predicted by droplet theory. 

For low temperatures the discreteness of the lattice becomes important. It has been 
shown |58H63| that for sufficiently low temperatures the lifetime of the metastable phase for 
the nearest-neighbor square-lattice Ising ferromagnet with Hamiltonian given by Eq. ( |T0| ) is 
given by 

ln(L 2 r) = 8(3l c - 2(3\H\(l 2 c -l c + l), (33) 

where the size of the critical droplet is l c = \2/\H\], and the notation \x] denotes the 
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smallest integer greater than x. This result is restricted to 2/\H\ not being an integer and 
to \H\ < 4. For this result to be valid the temperature must be at least low enough and 
the lattice size large enough to insure that the system is in the SD region. Differentiating 
Eq. ( p3|) with respect to gives 

TA = 2H 2 (l 2 c - l c + 1) . (34) 



V. NUMERICAL RESULTS 

In this section we present extensive results of the mean-field dynamics introduced in 
Sec. m for the relaxation behavior of the metastable phase in the nearest-neighbor Ising 
ferromagnet on a square lattice with the Hamiltonian given in Eq. flTOl). Periodic boundary 
conditions are used throughout. 

First, we consider T=0.8T C (j3=0. 55086. . .) for direct comparison with recent MC results 
0. We obtained Fo(m) for the entire range of— 1 < m < +1 for system sizes up to 
64x64, using the method outlined in Sec. [IV] and utilizing the Ising spin-reversal symmetry 
F (m) = F (—m). For H ^ 0, F(m) was obtained from Eq. ([[]). Once F(m) was obtained, 
the Markov transition probability matrices Wj and W 2 for the two dynamics MFD1 and 
MFD2 were constructed through the procedure outlined in Sec. ||. We show the results 
from MFD1 first and discuss MFD2 later. As in Ref. || the initial state is chosen as m=+l 
with H<0, and an absorbing barrier is put at m cut =0. In Ref. || different values of m cut 
were also used. As long as m cu t is sufficiently far away from the metastable value of m that 
the largest droplet must already be supercritical, the precise value of m cut is not important 
for weak fields 0. 

Using the absorbing Markov chain method discussed in Appendix O, the average first- 



passage time (r) to m=m cut and the standard deviation cr T =J (r 2 ) — (r) 2 were obtained by 



matrix inversion using the subroutine tridag from Ref. ||64|| . Our values of r were divided 
by L 2 to give all times in units of MCSS. 
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In Fig. [3] we show the relaxation time, defined as the first-passage time to m cut =0, ob- 
tained from MFD1 for L=32 and 64. For comparison, MC results for the standard Metropolis 
algorithm with spin updates at randomly chosen sites for L=32 and 64 are also plotted. We 
find that MFD1 and the microscopic MC simulations give qualitatively similar results for 
the relaxation time. 

In Fig. |] we show the slope, A, of the data in Fig. || with respect to whose 
asymptotic values in the SD and MD regions are given in Eqs. ( 2§D and (|30|) , respectively. 



We clearly see four distinct relaxation regions for L=64. The SF region, \H\ > Hmfsp 
(-^mfsp(0.8T c ) ps 0.75 fl2fl), contains the sharp peak for very small in Fig. £|(b). 



(Fig. |K a) is analogous with the schematic Fig. |2], in which the various field regions and 
crossover fields are indicated.) The region 0.2<|if |<0.75 corresponds to the MD region. As 
| if | is lowered, the crossover from the MD region to the SD region is signalled by a sudden 
rise in A. In Fig. |](a) we define i? max and H m [ n as the fields at which A has a local maximum 
and minimum, respectively. As \H\ is lowered further, A plunges towards zero, signalling the 
CE region. From our numerical data for MFD1 we find this crossover field i?THSP a t about 
ps 40 for L=64 and [Hi' 1 ps 20 for L=32. (See Fig. |(b).) This is consistent with 
the relation .Hthsp oc L~ l given in Sec. |m[ For further discussion of the thermodynamic 



spinodal and its relation to the droplet-to-slab transition, see Sec. [VI. 

Using the numerically exact value f43| , |4"8|| of 5, we find our results consistent with 6+c~2 
for MFD1, based on the data in the SD region. The expected value for dynamics that can 
be described by a Fokker-Planck equation is b+c=3. The estimate 6+c~2 for MFD1 at 
0.8T C should be taken with extreme caution, since the asymptotic region for L=64 appears 
to be quite small. For a more reliable estimate, we would need results from larger systems 
at this temperature. Since we have difficulties in obtaining F(m) for larger system sizes, we 
instead tested Eq. ( |28|) at lower temperatures, assuming that b+c does not depend on the 
temperature. Although the accessible system sizes are smaller for lower temperatures, due 
to the difficulties in estimating S(E, M) for all values of E and M, this approach turns out 
to be a more reliable way with our dynamics to estimate b+c. These studies at lower T are 



also consistent with 5+c~2, as discussed below. 

Although the results from the mean-field dynamic are in good overall agreement with 
those from the MC simulations, the quantitative estimates of the intercept and the slope 
in Fig. |](a) are less satisfactory for the MD region. Since the asymptotic region, in which 
Eq. ( |5UD is valid for MC simulations at 0.8T C , has been shown to be relatively narrow, even 
for L=720 H, results for much larger system sizes than L=64 would be necessary to provide 
a satisfactory test of Eq. (|30| ) for the mean-field dynamics at this temperature. At lower 
temperatures we are nevertheless able to estimate below that 6+c~2 for MFD1 in the SD 
region, and this gives the slopes drawn in Fig. |](a) for both regions. In contrast, it was 
demonstrated in Ref. for systems with L<720 that Metropolis MC with spin updates at 
randomly selected sites gives 6+c~3, in agreement with theoretical expectations P JT0| , p7| . 
The apparent consistency with 6+c~2 of the Metropolis MC results for L=64, which are 
also shown in Fig. |](a), is therefore clearly due to finite-size effects. 

In Fig. [| we show the field dependence of the slope A obtained from MFD1 and MFD2 at 
T=0.8T C . Both are in qualitative agreement with Fig. ^|. The SF region, where A decreases 
to zero for large \H\, is not shown. Figure ^ provides some insight about the relevance 
of F(m) and the possible artificial results from the mean-field dynamics. The fact that 
Fig. H is in qualitative agreement with Fig. ^| indicates the importance of F(m) for the 
dynamics. However, as one might expect from the sensitive dependence of the kinetic- 
prefactor exponent c on the details of the dynamic, which was demonstrated in Ref. ||, values 
of b+c obtained from mean-field dynamics do not correctly reflect that of the underlying 
microscopic dynamic. In the SD region, A approaches the exact value /3S with a slope that 
depends on the details of the particular dynamic. For MFD1 and MFD2, b+c is about 
2.0 and 3.7, respectively. In the MD region for both MFD1 and MFD2, we find only 
qualitative agreement with droplet-theory predictions. As will be explained in Sec. [VT], it 
is expected that only in the SD region can quantities related to the metastable phase be 
reliably extracted from F(m). In the MD region A depends strongly on the details of the 
particular dynamic, and the dynamics associated with the droplet growth is not taken into 
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account correctly in F(m). Note that the MD region is fairly easily accessible by standard 
MC simulations, since r there is rather small. 

In Fig. ^ we show the field dependence of the relative standard deviation of the lifetime, 

r = a T /(r) , (35) 

which provides additional information about how the metastable phase decays. If the decay 
of the metastability involves a single Possion process of forming one critical droplet, as is the 
case in the SD region, we expect r~l. In the MD region, on the other hand, one needs to 
consider many independent Poisson processes. By partitioning the system into [L/Rq) 2 ~^>1 
cells of volume proportional to Rq, one gets r oc Rq/L, where R is given by Eq. fl3T|) |5|,|65 



A more rigorous argument, based on the two-point correlation function ||57|| , can be found 



in the appendix of Ref . . An estimate for the crossover field -Hdsp between the SD and 



MD regions can be chosen as the field i/1/2, for which r=l/2 [|]||. For further discussion of 
the dynamic spinodal in mean-field dynamics, see Sec. [VT]. 

In Fig. 0(a) we show the temperature dependence of Hi/ 2 for L=24. The MC values lie 
between the estimates from MFD1 and MFD2. Again we considered the standard Metropolis 
dynamic with spin updates at randomly chosen sites. The MC simulations in this case 
were accelerated by using the method of absorbing Markov chains. This new MC method 
generalizes the n-fold way algorithm [0,^] and gives large CPU-time savings for low 
temperatures without changing the underlying dynamics. 

Analytic estimates for -£Tdsp & t l° w temperatures can be obtained as follows. For suffi- 
ciently low temperatures there exists a field 2<\H\<A such that a single overturned spin is 
a supercritical droplet |58H63|. We define T\ as the average time before a single overturned 



spin appears. Further, we define r 2 as the average first-passage time from the state with a 
single overturned spin to the absorbing state with magnetization m cut =0. For 2<|H|< 4 and 
low temperatures, the processes that determine r 2 are deterministic, and for both MFD1 and 
MFD2, t 2 is of order unity. Using the free-energy difference between the state with a sin- 
gle overturned spin and the metastable state with no overturned spins, 2(3\H\ — 8/3 + In A", 
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one can obtain t\ for MFD1 and MFD2 from Eqs. (§) and (PD , respectively. We get 
ri oc exp(8/3 - 2/3\H\ - InN) for MFD1 and n oc exp(8/3 - 2^[i?| - 21nJV) for MFD2. 
Note that we here give t\ in units of MCSS and that N=L 2 . The waiting time n can be 
either large compared to T2, corresponding to stochastic decay, or small, corresponding to 
deterministic decay. Therefore, the estimate for i^DSP can be obtained by setting t x to be 
of the same order of magnitude as r 2 , which leads to 

4-# DSP = T(lnL + Cl ) for MFD1 (36a) 
4 - H DSP = T(2 In L + c 2 ) for MFD2 , (36b) 

where C\ and c 2 are non-universal constants. For MC simulations, T\ oc exp(8/3 — 2(5\H\ — 
In TV) as for MFD1. Simulations in which the update sites are chosen sequentially give r 2 
on the order of unity, so that i^DSP is given by Eq. (|36a|) for these microscopic dynamics as 
well. When updates are performed at randomly chosen sites, however, r 2 increases with L. 
The radial growth velocity v is independent of L and H for 2<\H\< 4 and T<Cl, so that 
vr 2 = O(L), and t\ = r 2 then gives 

4 - # DSP = T (J In L + c 3 ) . (36c) 

In Fig. we compare these analytic low-temperature estimates for i^DSP with i?i/ 2 as ob- 
tained both from our mean-field dynamics and from microscopic MC simulations. Excellent 
agreement for the linear T dependence is demonstrated in Fig. 0(a) and for the logarithmic 
L dependence in Fig. 0(b). 

In Fig. |8| we show the field dependence of the slope A for L=24 as a function of \H\ for 
T=1.0, 1.1 and 1.2 in the SD region. (Note that T c =2.269- • •.) We observe that, as expected, 
the size of the asymptotic SD region increases significantly as T is lowered. The size of the 
MD subregion, which is not shown, is very compressed for these low temperatures, since the 
entropic factor is much less important in F(m). From Fig. |8| we confirm that 6+c=2.0(l) for 
MFD1 and 6+c=3.7(2) for MFD2. To obtain these estimates of b+c and /?H, we performed 
least-square fits to Eq. (|30|) considering only data points in the asymptotic SD region in 
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the following way. We first discarded data in the CE region. We then estimated the slopes 
and the intercepts as we successively removed data points from the strong-field end of the 
\H\ interval. Since the asymptotic region for Eq. ( p0|) corresponds to i^THSP< 1-^1^1; our 
estimates for the slopes and the intercepts approached constant values with small fluctuations 
as we removed data points. Due to the uncertainty in S(E, M) it is not easy to obtain a 
systematic error analysis. We therefore estimated the errors from the fluctuations of the 
estimates in the asymptotic region. The resulting estimates for /3H from MFD1 are slightly 
different from those of MFD2, as is shown in Fig. [8[ The final error bars in our estimates 
for (3E include this effect. Our estimates for (3E, 4.95(2), 4.00(1), and 3.20(2) for T=1.0, 
1.1, and 1.2, respectively, are within 1% of the exact values [pjPI 4.942, 4.004, and 3.217, 
respectively. 

Since the nucleation rate for a single droplet is independent of the system size, the 
average metastable lifetime in the SD region should be proportional to L~ 2 , as indicated by 
Eq. (P?D. We checked this result at T=1.0 and \H\~ 1 =2A7 by fitting the lifetime to the form 
(r)ocL~ Q for L=12, 14, 16, 20, and 24 (in the SD region for all the values of L used). From 
this we obtained a=2.08(3), in reasonable agreement with the theoretical result. At higher 
temperatures the agreement is less convincing, so that at 0.8T C we found a~3 for the same 
values of L. However, we believe this is a finite-size effect which becomes more pronounced 
at higher temperatures. 

It is also possible to construct mean-field dynamics that interpolate smoothly between 
MFD1 and MFD2. In analogy to Eqs. (§) and (pTTj) the transition probabilities of such a 
dynamic can be written as 
/ n \ 7 

W 7 (n,n+1) = A\1-—J exp[min{0,(F(n) - F(n + 1) + 7 (£(n) -S(n + 1)))}} (37a) 
W 7 (n, n - 1) = A C^Y exp [min {0, (F(n) - F(n - 1) + 7 (S(n) - S(n - 1)))}] (37b) 
W y (n, n) = 1 - W 7 (n, n + 1) - W 7 {n, n - 1) . (37c) 

The positive constant A only needs to fulfill the requirement that W 7 (n, n) >0 for all n, and 
it is otherwise unimportant since it only redefines the overall timescale of the process. In 
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principle, the only restriction on 7 should be 7>0. For 7=1 and A=l the process reduces 
to MFD1, while for 7=0 and A— 1/2 it reduces to MFD2. Since both the prefactors b+c in 
the SD region and the L dependence for the low-temperature strong-field behavior of -ffosp 
are different for MFD1 and MFD2, it is expected that they both change continuously with 
7. Thus it may be possible to tailor a mean-field dynamic to give the desired prefactor and 
Hdsp dependence to match a particular microscopic dynamic. For MC with sites selected 
sequentially, MFD1 has the desired values of these two quantities. If 7=1/2, then the T 
dependence of i^osp will be given by Eq. (|36cj ), which agrees with MC with randomly selected 
sites. The mean-field dynamic with 7=1/2 might also be expected to give a prefactor closer 
to the value for MC with randomly selected sites, b+c=3 0, since this lies between the 
values we have obtained for MFD1 and MFD2. 

In Fig. |9| we show TA for the very low temperatures T=0.4, 0.2, and 0.1. Discrete-droplet 
results for low T |58|- |63|1 are shown as a set of parabolic arcs given by Eq. ([34]). For T=0.4 
we observe clear oscillatory behavior with \H\, even though quantitative agreement with 
the discrete-droplet limit has not set in yet. For T=0.2 and 0.1 we observe increasingly 
good quantitative agreement between the MFD1 results and the discrete-droplet limit. For 
T=0.4, (r) is on the order of 10 7 and 10 26 for |if|=1.0 and |i?|=0.3, respectively, as obtained 
from MFD1 with m cut =0. The corresponding numbers for T=0.2 (0.1) are 10 18 (10 39 ) and 
10 59 (10 123 ). 

MC results by one of us provide corroboration for the MFD1 results. This oscilla- 
tory behavior with \H\ is due to the corrugation of the free energy near the saddle point, 
illustrated in Fig. 0(b). This corrugation is due to the discreteness of the lattice. 

VI. RELATIONS TO EARLIER WORK 

The relationships of the detailed shape of F(m) to two-phase equilibria and nucleation 
barriers in finite systems have previously been studied by Binder and coworkers fl8|-|2"0"fl . 
Although their work covers a wide range of temperatures below T c , the bulk of their analysis 
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is concerned with the critical region, in which the correlation lengths become comparable to 



the other length scales, L, R , and R c (see also Ref. ]6"Bf| ). Since F(m) can be viewed as a 
potential for the stochastic processes that define our macroscopic dynamics, it might have 
been expected that those studies also would have revealed two L-dependent spinodals with 
the same scaling behaviors as .£/thsp an d H DSP . However, they instead reported a single 
spinodal at a field proportional to L^ 1 [19|. 



To clarify the relation between our results and those of these earlier studies, we follow 
Furukawa and Binder by calculating 



-d\-i dF o( m ) 



h(m) = (/^T 1 ^^ • (38) 

This quantity has the dimension of a field and is the equivalent, for an Ising model in the 
fixed-m ensemble, of the expectation value of the chemical potential for a lattice gas in the 
fixed-density ensemble [19|. In Fig. [1^ we show h(m) for L=24, 32, 64, and 96 at 0.8T C [69 



(For L=96 the simulation was only performed for |m|>0.63.) From Eqs. ([]]) and ( |38"D it is 
seen that 

= (3L d (h(m) - H) . (39) 

am 

Thus the extrema of the field- dependent F(m) occur where the applied field H equals h(m). 
For a negative H between zero and the minimum of h(m), there are three such extrema: the 
stable minimum at m~— m eq (T), the metastable minimum at m^+m eq (T), and an unstable 
maximum at some m between and m sp , the magnetization corresponding to the minimum 
of h(m). If dynamically relevant information for nonzero H is to be deducible from the 
functional form of F(m), this must mean that the spatial configuration characteristic of the 
nonequilibrium saddle point of the relaxing system is well approximated by the equilibrium 
configuration corresponding to the same value of m in the fixed-m ensemble. It is not 
unreasonable to expect that this should hold, at least for sufficiently large systems that the 
corresponding applied H is weak. 

Following the reasoning outlined above, we start by considering the thermodynamic 
spinodal. The droplet-theoretical arguments in Sec. [IV] indicate that the THSP should 
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be located in the field range where R c becomes on the order of L, so that the droplet 
must compete with slab configurations about being the "true" saddle point. It has been 
shown by Leung and Zia ||2] that for Ising models with periodic boundary conditions in the 
fixed-m ensemble, a first-order phase transition between the equilibrium droplet and slab 
configurations occurs at a temperature dependent magnetization, m c (T). For \m\<m c the 
equilibrium configuration is a slab. The volume fraction corresponding to m c (T) in the limit 
L-^oo, which was obtained in Ref. [ 42| , can be written as |49| 



The vertical, dashed lines in Fig. |10| mark m c (0.8T c ). If one ignores entropy effects due 
to the center-of-mass positions of the droplet and slab, as well as capillary waves on their 
surfaces (which are responsible for the \H\ b power-law prefactor in the nucleation rate 
the finite-size displacement of m c is proportional to L v WA - The vertical arrows in Fig. 
indicate the corresponding values of m c for L=24, 32, and 64, as obtained from Fig. 3 of 
Ref. These values monotonically approach the infinite- L limit for m c . From Fig. |I"0|(a) 



it appears that the shape of h{m) evolves towards a step discontinuity as L increases. From 
the data collapse in the plots of Lh(m) in Fig. |10](b), it is seen that the magnitude of this 
discontinuity is proportional to L~ x and so corresponds to a finite discontinuity in a first 
derivative of the free energy per unit interface area, L~( d ~^ -Fo(m) , i n the fixed-m ensemble. 
This is consistent with the identification of the droplet-to-slab transformation as a first-order 
phase transition f4"2fl . 

Using 0c (T) from Eq. (|I]) in Eq. fl2"5| ) for if TH sp(0), we obtain an estimate for H THS p 
as the field at which the free energies of the critical droplet and a system-spanning slab are 
degenerate: 

1 ( d SfT)\^ 

ff T H S P = ffx„sp(« = I (^j ■ (41) 

The horizontal, dashed line in Fig. |T0|(b) represents this estimate of L.£/thsp- It corresponds 
excellently to the magnitude of the incipient step discontinuity in Lh(m). The L dependence 
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of Hthsp at 0.8Tc is shown in Fig. [11], together with the results of MFD1 for L=32 and 64. 
The MFD1 points were estimated as the smallest fields for which A(T, H)=f3'E(T) (see 
Fig. |J). The agreement is good, and we believe the small discrepancy is a finite-size effect. 

The narrow minimum in h(m), which occurs at an L-dependent magnetization m sp (L), 
slightly smaller than the equilibrium magnetization m eq , signifies the inflection point in 



Fo(m). As pointed out by Furukawa and Binder |19|], this magnetization also corresponds 
to a change in the equilibrium configuration: closer to m eq the system is uniform, whereas 
closer to m c a single droplet of the opposite magnetization is precipitated. The volume 
fraction occupied by this droplet is given by the lever rule. As seen from Fig. [13(b), the 
minimum value of h(m), h(m sp ), does not vanish as L^ 1 . 

It is tempting to identify the disappearance of the single-droplet saddle point, which 
occurs at h(m sp ), with the dynamic spinodal, H DSP . For all the system sizes studied here, 
we find that \h(m sp )\ indeed lies close to other estimates for i^osp, such as Hyz obtained 
both from MFD1 and from microscopic MC simulations. These comparisons are illustrated 
in Fig. |TT[ Obviously, the range of L used in the present study is too narrow to obtain the 
scaling relation with any degree of certainty. However, the values of \h(m sp )\ at different L 
also agree well with the estimate for i^DSP obtained in Ref. |5| by fitting the proportionality 
constant in the relation LocRq with Rq given by Eq. (|31~1) to MC data for L between 64 
and 720. We find this noteworthy, considering that this analytical expression for H DSP 
is based on explicitly dynamical arguments, whereas \h(m sp )\ is obtained from a strictly 
equilibrium calculation. 

To account for this observation, we suggest that m sp corresponds to the volume fraction 
at which entropy effects make the free energy of a configuration consisting of two droplets 
lower than that of the single-droplet configuration [[0J. We assume that the single droplet 
is replaced by two identical droplets, each with half the volume of the original droplet. We 
neglect corrections to the droplet free energy, including those arising from surface excitations, 
which correspond to the power-law pref actor \H\ b in Eq. (|2"6|). We only consider the entropy 



contributions due to the droplets' center-of-mass positions, and we neglect excluded-volume 
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effects. This gives the approximate zero-field free energies 

d-l 

F o (1) (0) « F (m cq (T)) + L d ~ l df3 (^^) ' H(T)3 - dlnL (42a) 
for the single-droplet configuration and 



d-l 



F o (2) (0) « F (m eq (T)) + ifrL^W " S(T)* - 2dlnL (42b) 

for the two-droplet configuration. Equating an d we find 



d-l / TlnL 



d 
d-l 



SP 2m eq (T)5(T)^L d V2^-l, 



(43) 



Inserting sp into the single- droplet approximation for obtained from through 

Eq. (PS|), which is identical to the expression for if TH sp(</>) given in Eq. (^), we obtain 

/(2^-l)/5S(r)\^ 



Comparing this result with Eq. ( p2|) for i^DSP we n °t e that, except for a ci-dependent numer- 
ical constant of order unity, it has the same asymptotic dependences on L and T as -£/dsp- 
In particular, |/t(0 sp )| does not vanish as but rather much more slowly as (lnL)~~. 
As noted previously and illustrated in Fig. |TT], corrections to this asymptotic behavior for 
Ht>spi due to the power-law prefactors in the nucleation rate, are very substantial. Since 
analogous corrections were ignored in the approximate derivation of \h(<j) sp )\ given here, we 
expect similarly large corrections to apply to it. This is in agreement with the numerical 
results shown in Figs. [10] and [11]. We further note that the precise value of the numerical 



coefficient in our approximate expression for the asymptotic value of \h((f) sp ) \ is the result of 
our choice to consider only a separation of the single droplet into two equal droplets. A more 
careful calculation ought therefore to give a different coefficient. However, the factor (d+1) 
in the denominator of i^DSP results specifically from the simultaneous nucleation and growth 
processes which give rise to the expression for the lifetime in the MD region, Eq. fl2"9"|). It is 
therefore unlikely that further improvement of the equilibrium calculation should yield the 
same factor in \h(<j) sp )\. 
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For |iy|>|/i(m sp )| the bulk free energy displays no saddle point, but its derivative with 
respect to m has a minimum at m sp . As a consequence, mean-field dynamics should show a 
slowing-down near m sp , as observed in the present study. However, F(m) does not contain 
information about the complicated multi-droplet configurations that dominate the dynamics 
in this region of relatively strong fields. Consequently we expect the agreement between 
microscopic and mean-field dynamics to be only qualitative in the MD region. 

VII. DISCUSSION 

In this paper we have introduced a class of macroscopic mean-field dyamics, and have 
studied in detail two members of this class which we call MFD1 and MFD2. We have 
demonstrated that these macroscopic mean-field dynamics replicate many of the qualitative 
and quantitative features of the relaxation behavior of the metastable phase of the two- 
dimensional nearest-neighbor Ising model. As a function of the external field H, four distinct 
regions of relaxation were observed in agreement with recent microscopic MC studies . In 
the single-droplet region at low temperatures, the leading exponential term in the relaxation 
time (r) was obtained to within 1% of the exact value. We also obtained temperature 
independent estimates for the prefactor exponents b+c for the two mean-field dynamics. At 
very low temperatures we observed an oscillatory behavior in A(if)=dln(r)/d|iJ| _1 with 
respect to \H\, in agreement with discrete-droplet theory. In the low-temperature limit, 
simple theoretical estimates of the dynamic spinodal field i^DSP; i n terms of the temperature 
and the system size, were obtained for various dynamics. Our numerical studies provide 
excellent agreement with these predictions. 

The mean- field dynamics are constructed with only the following conditions: (1) locality 
in the value of the order parameter m, and (2) the correct equilibrium distribution obtained 
from the order parameter of the microscopic model. These two conditions constitute the 
minimum requirements for any local dynamic. The reasonable results obtained from the 
macroscopic mean-field dynamics may be somewhat surprising in view of the fact that no 
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microscopic information is directly included. However, we believe the relative success of our 
dynamics becomes more understandable if one views them as approximations to the macro- 
scopic dynamic which can formally be constructed by using a projection-operator formalism 
to project the microscopic dynamic onto a master equation for the one-dimensional order 
parameter distribution, as discussed in Appendix The dynamically relevant information 
which is retained in F(m) following this projection correctly describes the droplet configu- 
rations that provide the rate-determining steps in the decay process. In particular, these are 
the droplet and slab configurations that are important near the thermodynamic spinodal, 
the single-droplet configuration characteristic of the SD region, and the breakdown of the 
single-droplet configuration into a uniform "gas" of microscopic fluctuations that takes place 
near the dynamic spinodal. These aspects were discussed in light of earlier work in Sec. [V|. 
As a consequence, the mean-field dynamics produce excellent numerical estimates for both 
the thermodynamic and the dynamic spinodal fields. 

Comparison between the two proposed dynamics, MFD1 and MFD2, provides some 
insight into the dependence of the lifetime r on the detailed dynamics. Most of the char- 
acteristic behavior of A predicted by both continuous- and discrete-droplet theory |5U5g|] is 
expected to hold for different local dynamics. However the influence of the detailed dynamic 
is reflected in the prefactor exponent b+c. The difference between A (if) as obtained from 
MFD1 and MFD2 in the deterministic region suggests that A (if) in this region depends 
strongly on the particular dynamic. 

Generally, when the relaxation time is long, the system spends more time exploring phase 
space, and therefore the dynamic is more strongly subject to the restricted bulk free energy 
F(m), as indicated in this study. When the relaxation time is short, the details of the 
particular dynamic are more important than the bulk free energy. From the field-theoretical 
point of view, when the most probable trajectory from the metastable phase to the unstable 
saddle point is sharply defined, the relaxation time is usually large, and regardless of the 
details of the particular dynamic, one needs to consider only the trajectories near the most 
probable one in order to study the decay of the metastable phase. When the probability 
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distribution over trajectories is not very sharp, the lifetimes are usually short and the details 
of the particular dynamic play important roles. 

The biggest advantage of macroscopic mean-field dynamics is that they can provide data 
for low temperatures and weak fields (in the single-droplet region), for which the average 
lifetime of the metastable phase is too long to be measured with standard MC algorithms. 
They also provide data for arbitrary temperatures and fields, allowing one to obtain accurate 
estimates of derivatives. However, the system sizes for which mean-field dynamics can be 
applied are rather limited by computational constraints. 

In summary, we have presented a method to study the relevance of the equilibrium 
properties of a model to the dynamical relaxation of metastable phases. The macroscopic 
dynamics are designed using only the minimal requirements of locality in the relevant order 
parameter and the correct equilibrium free energy projected on that order parameter. Ex- 
tensive applications to the two-dimensional nearest-neighbor Ising ferromagnet on a square 
lattice provides convincing evidence that the characteristic behavior of the dynamical re- 
laxation of the metastable phases is largely determined by the restricted bulk free energy 
F(m). We believe that our approach can benefit studies of relaxation phenomena for other 
systems as well. 
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APPENDIX A: RELATION OF OUR MACROSCOPIC MAGNETIZATION 
DYNAMICS TO THE MICROSCOPIC SPIN DYNAMIC 



In this Appendix we briefly consider the formal relationship between our macroscopic 
magnetization dynamics, MFD1 and MFD2, and the underlying microscopic dynamic repre- 
sented by the single-spin-flip Metropolis algorithm with updates at randomly selected sites. 
The formal framework for the discussion is the Nakajima-Zwanzig |^,^] projection-operator 
formalism for the master equation, which is equivalent to Mori's |27| projection-operator for- 
malism for the equations of motion of observables P8|J23[]. We adapt the standard discussion 
(see, e.g., Refs. p9| , |30[| ), which considers a deterministic microscopic dynamic governed by a 



quantum-mechanical or classical Hamiltonian, to the case where the microscopic dynamic is 
a discrete-time Markov process (possibly derived from a deterministic dynamic at an even 
more microscopic level). 

For an iV-site kinetic Ising model, the microscopic probability density at time k is a 
2 N - dimensional column vector p(k) which evolves in time according to the equation 

p{k+l) = Wp{k) , (Al) 

where W is the matrix of microscopic transition probabilities. The probability distribution 
over the "relevant" macroscopic variables at time k, X(k), is obtained from p(k) through 
the action of a projection operator V, 

X(k) = Vp{k) . (A2) 

Although the dimension of X(k) is 2^, if the only macroscopic variable considered is the 
magnetization, the dimension of the "relevant" space in which X(k) has nonzero components 
is N+l. By using Eqs. (|A1|) and (|A2|), one can write the equation of motion for X{k) as 



X(k+1) = VWp(k) = VW(V+Q)p(k) 

= VWX(k)+VW[Qp(k)} , (A3) 
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where Q=l— V is the projection operator onto the "irrelevant" orthogonal complement to 
the relevant space. The first term in the second line of Eq. ( |A3| ) corresponds to a Markov 
process for X(k), whereas the second term contains non-Markov contributions which can be 
formally evaluated as follows. 

Operating on Eq. ( |A1| ) with Q one obtains 

Qp(k) = QWp(k-l) = QW(V+Q)p(k-l) 

= QWX(fc-l) + QW [Qp\k-1)\ , (A4) 



which is inserted in Eq. (|A3|). Iterating this procedure a total of k times, one obtains the 
final result: 

k 

X(k+1) =VWX(k) +PW^[QW]'X(fc-I) + VW[QW] h Qp(0) . (A5) 

1=1 

The first term on the right-hand side corresponds to a Markov process in the relevant sub- 
space. The non-Markovian second and third terms represent memory about the relevant 
variables at earlier times, propagated through the irrelevant subspace, and specific informa- 
tion about the initial state of the irrelevant variables, respectively p0|| . The third term can 
usually be ignored, at least after a short initial period. 

In standard applications of projection-operator techniques, the quality of the resulting 
approximation depends on the choice of the relevant macroscopic variables. The approach is 
most useful whenever there is a large separation between "fast" and "slow" timescales, and 
one usually attempts to include all the slow variables in the relevant subspace. This ensures 
that only variables with short correlation times contribute to the memory effects, which can 
then often be ignored or approximated by a rapidly decaying function. The (most obvious) 
slow variables are determined by macroscopic conservation laws or by spontaneously broken 
symmetries in the corresponding isolated system [2^ . 



In the present work we have considered the stochastic time evolution of the magneti- 
zation (our relevant macroscopic variable) as a Markov process defined by the transition 
probability matrices Wi (for MFD1) or W2 (for MFD2). In doing so we have performed 
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a Markov approximation equivalent to ignoring the memory effects and using Wi and W2 
as approximations for the matrix obtained by contracting VW in Eq. (|X5|), so that its di- 
mension becomes N + 1. The macroscopic slowness of the magnetization is related to the 
spontaneously broken symmetry between the two ferromagnetic phases for H=0 below T c . 
By virtue of energy conservation in the corresponding closed system, the other obvious slow 
macroscopic variable is the total energy. In relegating it to the irrelevant subspace, mainly 
for computational convenience, we have most likely ignored non-negligible memory effects. 
By considering only a single relevant variable we also have excluded nonlinear interactions 



between macroscopic variables 2S 



The satisfactory agreement between our approximate macroscopic Markovian dynamics 
and the MC simulations of the full microscopic dynamic indicates that the approximations 
made in the present work are quite reasonable. Nevertheless, the above discussion indicates 
that by including the total energy as a second relevant macroscopic variable in our Markovian 
mean-field dynamics, we could reduce the importance of the neglected memory effects and 
allow for nonlinear interactions between relevant variables. We believe further significant 
improvement of the agreement between the approximate, macroscopic dynamics and the 
underlying microscopic dynamic could be achieved in this way. 
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APPENDIX B: ABSORBING MARKOV CHAINS 



In this Appendix we briefly discuss the application of absorbing Markov chains |32j to 
obtain the expectation value and the variance of the first-passage time for escape from 
the metastable phase. For a given transition probability matrix W(n,n'), simple matrix 
calculations provide the expectation value and the standard deviation of the first-passage 
time from one state to another. 

Let us consider a random walker which moves between states < n, n' < N with tran- 
sition probabilities W(n,n'). Starting from an arbitrary initial state X(0), the probability 
density after k time steps, X(k), is 

X{k) = W k X(0) . (Bl) 

Without absorbing states, the probability of a Markov chain is conserved, that is, 

J2W(n,n') = l (B2) 

n' 

for all n. After k time steps the walker is still in some state, that is 

e r W k X(0) = 1 , (B3) 

where e 1 ' = (1, - • • , 1, • • • , 1) and the superscript T denotes the transpose. 

Now we place absorbing states at i > n cut . Once an absorbing state is reached, the 
walker is absorbed and the Markov chain terminates. Let T be the n cu t x n cnt submatrix of 
W that contains the transition probabilities between the n CVLt transient states. The analog 
of Eq. ( B2|) is not satisfied for T. The probability that the walker is absorbed at time k is 
e r (T k ~ 1 — T fc )X(0). The average first-passage time to the absorbing states is 

oo 

(r) = e r k(T k - 1 - T k )X(0) = e T NX(0) . (B4) 

fc=i 

The fundamental matrix is defined by N = (I — T) _1 , where I is the identity matrix. 
Similarly, the second moment of the first-passage time can be obtained from |32j 
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(r 2 ) = e*(2N 2 - N)X(O) . (B5) 

In this work, r was divided by the total number of sites N, so that all times are given 
in units of MCSS. 
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FIGURES 

FIG. 1. (a) The zero-field bulk restricted free energy, Fo(m), of the nearest-neighbor Ising 
ferromagnet on a square lattice at T=0.8T C =1.8153. . . for L=24. (b) Same as (a), but at 
T=0.2~0.0881T C . Notice the difference in the scales along the y-axis in (a) and (b). The saw- 
tooth-like behavior of Fq(iji) at T=0.2 is due to the discreteness of the lattice. The vertical arrows 
marked m c indicate the exactly known magnetization for which the most likely configuration with 
the given magnetization changes from a slab (for \m\ < m c ) to a droplet (for \m\ > m c ) in an 
infinite system | |42fl . 

FIG. 2. Schematic plot of A.(H), defined in Eq. (pl|), for a two-dimensional Ising ferromagnet. 
The dynamic spinodal field -ffosp separates the stochastic and the deterministic regions. In the 
stochastic region (\H\<Hjysp), the relaxation time (r) is determined by the formation of a single 
critical droplet. Depending on the size of the critical droplet relative to the system size, the 
stochastic region is divided into the single-droplet (SD) and coexistence (CE) subregions. In the 
SD region the size of the critical droplet is smaller than the system size. The CE region is also 
characterized by a single nucleating droplet, but of a size comparable to the system size. The 
thermodynamic spinodal field Htusp separates the SD and CE regions. The deterministic region 
(\H\>Hdsp) is comprised of the multi-droplet (MD) and strong-field (SF) regions. These two 
regions are separated by the mean-field spinodal (MFSP). Droplet theory || predicts that the 
intercepts of the two straight lines are /3S and /3S/3. Their slopes are related as shown in the 
figure if one assumes that the radial growth velocity of the supercritical droplets is proportional to 
H. 

FIG. 3. The field dependence of the average metastable lifetime (r) for a two-dimensional 
Ising model at T=0.8T C . The lifetime (r) is estimated as the average first-passage time to m cut =0 
from the starting configuration m=+l with H<0. The lines correspond to the MFD1 dynamic for 
L=32 (dashed) and 64 (solid). The results of Metropolis MC simulations with random site updates 
for L=32 and 64 are marked by + and x, with 10 3 escapes from the m=+l state the statistical 
errors smaller than the symbol size. 
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FIG. 4. (a) The slope A(H), obtained from the data in Fig. [||, shown as a function of \H\. 
The thick solid curve is from MFD1 for L=64, and the symbols are from the standard Metropolis 
MC simulations with random site updates for L=64. The two straight lines are drawn with slopes 
b+c=2 and (6+c+2)/3= 4/3, using the exact value (43|,[4§ of /3H(0.8T C ) = 0.5062. . .. The horizontal 
arrows mark the exact values of /?E(0.8T C ) and /3S(0.8T c )/3. The asymptotic SD subregion seems 
to be very small for L=64. The SF subregion (\H\>2), where A(H) decreases to zero is not shown. 
We define H max and H m [ n as the fields at which A has a local maximum and minimum, respectively, 
(b) The slope A-(H), shown as a function of The two horizontal lines correspond to the 

exact values of /3S and /3S/3. The interpretations of the other lines and symbols are the same as 
in (a). 

FIG. 5. The field dependence at T=0.8T C of A(H) for the two different mean-field dynamics, 
MFD1 and MFD2. Qualitatively similar behavior is observed for both dynamics. Using the exact 
value of /3H, we estimate in the SD region 6+c~2.0 and 3.7 for MFD1 and MFD2, respectively. In 
the MD region, the lifetimes for MFD2 are shorter than for MFD1. Even though the system size is 
too small to discuss the asymptotic behavior in the MD region, the results from MFD1 and MFD2 
are not in good agreement with droplet theory in this region. 

FIG. 6. The field dependence of the relative standard deviation r=a T /(r) at T=0.8T C . The 
curves are from MFD1 for L=24, 32, and 64. The symbols are from the Metropolis MC simula- 
tions with random site updates for L=32 and 64. The MC data are from 3000 escapes from the 
metastable state near r~l/2, and at least 100 escapes from the metastable state for the other H 
values. 
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FIG. 7. (a) The temperature dependence of the estimate for the dynamical spinodal field 
-f^DSPi given by H\j2 for L=24. Note that the MC simulation results (Metropolis with updates at 
randomly selected sites, MR) fall between the MFD1 and MFD2 results. The estimates approach 
H=4 linearly in T with negative slopes, which depend on the system size L and a non- universal 
constant. The slopes of the straight lines are 5.7, 4.0, and 3.0. (b) The asymptotic slopes near 
the point (T, i? 1 / 2 )=(0, 4) as functions of L. We also show the MC results from the microscopic 
Metropolis algorithms with sequential (MS) and random (MR) spin updates. For comparison we 
show straight lines with the predicted slopes (2, 1.5, and 1 from top to bottom), obtained from the 
analytic low-temperature estimates for -ffDSP i n Eq. 



FIG. 8. The field dependence of the slope A(H) for L=24 and T=1.0, 1.1, and 1.2. The size 
of the asymptotic SD region increases as T is lowered. We estimate b+c = 2.0(1) and 3.7(2) for 
MFD1 and MFD2, respectively. The differences between the average values of /?E for MFD1 and 



MFD2 and the exact value of /?H [43,48 are less than 1%. Arrows indicate the exact values of /3E 



and the lines are from the fitted asymptotic behavior of A(H) in the SD region. 

FIG. 9. (a) The field dependence of TA(H) for MFD1 at T=0.4. The discrete-droplet 



result is shown as a set of parabolic arcs. The exact value of S [43],48| is indicated by the horizontal 
arrow. Simulation results for L=24 || from the Metropolis algorithm with spin updates at ran- 
domly chosen sites and 10 3 escapes from the initial state are also shown. The two vertical arrows 
indicate the estimates of -Hi/2 from MC simulation (left) and MFD1 (right), (b) Same as (a) at 
T=0.2. The results from the MC simulations and MFD1 agree quite well with the discrete-droplet 
|58|-[63]| result, (c) Same as (a) at T=0.1. The agreement between the three sets of results is close. 
Here the relaxation time (r) from MFD1 with m cu t=0 is on the order of 10 39 and 10 123 at |#|=1.0 
and 0.3, respectively. The MC estimate for (r) is on the order of 10 31 at |#|=1.35. 
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FIG. 10. The quantity h(m), defined in Eq. (p8|), shown vs magnetization m at T=0.8T C . 
The large x marks the spontaneous zero-field magnetization at this temperature. For detailed 



discussion, see Sec. VI. The vertical, dashed line in both panels marks m c , corresponding to the 
droplet-to-slab transition in the limit L— >oo, and the vertical arrows mark m c for the three smallest 
L studied (from left to right), (a) Shows h(m), highlighting the developing step discontinuity near 
m c , related to the thermodynamic spinodal, and the narrow minimum near m sp , related to the 
dynamic spinodal. From bottom to top the system sizes are L=24, 32, 64, and 96 (the latter 
only for m>0.63) (b) Shows L h(m), highlighting the L~ l scaling behavior of h(m) in the SD 
region between m c and m sp , as well as the much slower vanishing of the minimum value, h(m sp ). 
The horizontal, dashed line corresponds to L -£/thsp from Eq. (f4l|), and the dotted curve is the 
asymptotic single-droplet result for L h(m) obtained by applying Eq. (|38|) to Fq 1 ^ from Eq. fl42a| ). 

FIG. 11. "Spinodal phase diagram," showing the MD, SD, and CE regions in the (1/lnL, H) 
plane for the two-dimensional Ising model at 0.8T C . The lower, solid curve is -Hthsp from Eq. (fH|), 
and the data points (o) close to it are -Hthsp from MFD1 for L=32 and 64. The upper, solid 
curve is i^DSP from Ref . || , obtained from a one-parameter fit of the relation LocRq to MC data 
for L=64, 128, 256, 400, and 720. The short, dotted line indicates the asymptotic slope of -Hdsp- 
The data points represent different estimates for -Hdsp> obtained from MFD1, MC in this work, 
\h(m sp )\, and from MC for larger systems in Ref. ||, as indicated by the key in the figure. The 
quantities in the upper box are all estimates for -Hdsp- I n an cases, error bars are only given where 
the statistical error is larger than the symbol size. From right to left the data points correspond 



to L=24, 32, 64, 96, 128, 256, 400, and 720. See detailed discussion in Sec. VI. 
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